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We present a model of integrate and fire oscillators that move on a plane. The phase of the 
oscillators evolves linearly in time and when it reaches a threshold value they fire choosing their 
neighbors according to a certain interaction range. Depending on the velocity of the ballistic 
motion and the average number of neighbors each oscillator fires to, wc identify different regimes 
shown in a phase diagram. Wc characterize these regimes by means of novel parameters as the 
accumulated number of contacted neighbors. 
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1. Introduction 



Among the many emerging phenomena we observe in nature, synchronization is one of the most paradig- 
matic examples. It can be roughly understood as the collective dynamics of units whose internal 
state evolves p e riodic ally in time and when they interact tend to synchronize their internal variables 



Pikovskv et all |2001| ] . The achie vement of the final synchronized state (if a ny) strongly depends on the 



interaction pattern of the system Boccaletti et al. . 20061 ; Arenas et al\ . 20081 ]. Up to now, synchronization 
has been mainly analyzed in fixed topologies but we are witnessing the first evidences that l inks between 
agents can evolve in time (Olfati-Saber et al. . 2007 : Onnela et al . 2007 : Valencia et al. . 20081 ] . The case in 
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whic h this evolution of the network topology is an effect of the age nts mobility is a particularly interesting 



case 



Buscarino et al . 20061 ; Tanner et al. . 20031 ; Buhl et al . 2006 ]. The effect of this changing patterns of 



interaction on synchronization features has been anal yzed in different settings, for instance in chemotaxi s 
Tanaka , 20071 ] , mobile ad hoc networks Romerl . 2001 1, wire l ess se nsor networks Sivrikava fe Yener , 2004 ] , 



Uriu et aZ.I . [2010l ] . 



and the expression of segmentation clock genes 

In the recent literature, studies on synchronization in dynamically evolving complex networks have 
been mainly concentrated in th e case when the topology changes very fast. This i s the so-called fast 
switc hing approximation (FSA) (Belykh et al . 20041 : iFrasca et all . 20081 : Porfiri et al . 20061 : Stilwell et al. 



2006], which replaces the real interaction between agents by the "mean field assumption" that all agents 



interact with an effective strength that corresponds to the probability that any pair of agents are connected. 

Recently, it has been proposed a general fra mework of mobile oscilla tor networks where agents perform 
random walks in a two-dimensional (2D) plane [Fuji war a et oll . l2011bl ]. It has been shown that FSA fails 
when the time scale of local synchronization is shorter than the time scale of the topology change due to 
the agent motion. New behaviors arise due to the interplay between instantaneous network topology, agent 
motion, and interaction rules. This framework, that reduces to FSA when velocity is high enough, is valid 
for models whose evolution can be well app roximated by linear dynamics. This ac tually holds for models 
such as populations of Kuramoto oscillators [Kuramotol . ll984l : lAcebron et ~al\ . 120051 ] . whose evolution, after 
a short transient time, is very well describe d by a set of linear eq uations that can be solved in terms of 
spectral properties of the Laplacian matrix [Fujiwara et al . 2011al ]. 

In the present paper, we focus on a dynamical system, a population of integrate and fire oscillators 
(IFO), where linearization is not a good approximation, since the evolution takes place in two different time 
scales. One for the slow evolution of the internal state variables (the phase and the orientation) and the 
other for the fast interaction between the units (pulse coupling). During the last years it ha s been shown 
that t he interaction structure plays a fundamental role in the dynamics of IFO networks. IZillmer et al\ 



20091 ] observed different dynamical regimes due to network connectivity in a system formed by inhibitory 



integrate-and-fire neurons that were randomly connected. Also, the underlyi ng network structure can affect 
the speed with which the system reaches the synchronized state, as studied by Grabow et al. 201 Usually, 
IFO have been used to model neur al systems but we can also find some examples of applications in other 
fields, as for example in economy Erola et all 2011 ]. Models where the oscillators do not remain fixed, 
and the network of interactions changes with time can find a direct application in bio logical systems such 
as flashing fireflies, that in t ercha nge light signals while searching for potential mates [Mirollo &: Strogata . 
19901 : iRamirez-Avila et aZ.I . [20lTI ] . 



In the present case we will show that the interplay between agents motion and phase evolution towards 
a synchron ized state presents diff erent asymptotic behaviors, reminiscent of t he observation in Kuramoto 



oscill ators [Fuiiwara et all l2011bl ] and agents using communication protocols Baronchelli &l Diaz-Guileral . 
201 We identify, furthermore, the possible mechanisms in the different regions of the parameter space. 



The organization of the paper is as follows. In the next section we introduce the model. Then we show 
the results for different regions of the parameter space, velocity of the agents and range of interaction, 
and later we identify the different microscopic mechanisms that lead the system to a globally synchronized 
state. 



2. The model 

We p ropose a setting in which a population of N integrate and fire oscillators (IFOS) [Mirollo &: Strogatz 



1990] move at a constant velocity V in a bidimensional plane of size L with periodic boundary conditions 



Each agent has two degrees of freedom corresponding to an internal phase 4> € [0, 1] and orientation 
9i € [0, 2tt], both randomly set in an uniform manner at the initial configuration. 

The evolution of the system takes place on two different timescales. The slow timescale sets the pace 
at which the phases of the agents increase uniformly with period r, 

dt t 1 ' 

until they reach a maximum value of 1, when a firing event occurs. Then the phase is reset and the 
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oscillator is randomly reoriented. Upon this event at time t, the firing oscillator influences its nearest 
neighbors (oscillators at minimal distance) producing an update in their phases by a factor e: 

( Mt + ) = o 

= 1 =M <Pnn(t + ) = (1 + e)0 nn (i-) . (2) 
[0i(t+) G [0,2tt] 

The phase and orientation resetting corresponds then to the fast time scale. If the neighbor's phase update 
overcomes the phase maximum, another firing event is triggered and this process goes on repeatedly until 
all shots have ceased. At this point, the time t runs again until the next firing event occurs. The system 
is synchronized when we encounter in the system a succession of consecutive firing events (avalanche) 
equal to the system size N, since after this fact all oscillators will remain synchronized forever because 
all of them will have the same period r with or without interactions. For the sake of clarity we define the 
(discrete) time T, as the number of times a given oscillator (that we will identify with oscillator 1 in our 
computer simulations) has fired. This allows us to define T sync as the minimum number of (integer) cycles 
this reference oscillator takes to enter the synchronized state (i.e. the number of updates needed for an 
avalanche of size N to occur). 




Fig. 1. The model of interaction between oscillators, based on geometrical constraints. Only those within a distance R are 
affected by the firing of the central one. 



We propose a geometric condition for neighbor selection upon a firing event as shown in figure [TJ Every 
agent scans a circular area of radius R around it and shots the neighbors therein. We introduce a parameter 
r £ [0, 1] that indicates the fraction of the system available for interaction and relates both R, L variables 
and the average outgoing degree of the nodes of our evolving network 

tt/? 2 

(k out ) = (N- l)r. (3) 



L' 2 



0.02 ~ 0(1/N), 



Throughout this paper, we have used fixed parameters L = 100, r = 1, N = 50 and e 
while analyzing the explicit dependence on the mobility parameters, r and V . 

Before proceeding to show the results of our simulations, we need to note that the type of proposed 
interaction range in this sy s tem ha s been reported to show statistical properties similar to a continuous 
percolation [Fujiwara et al. . 2011b|, that i n the case of static oscilla tors occurs for approximately r c ~ 
4.51/(iV - 1) = 0.09 [Pall & Christensenl . 120021 : iBalister et all . 120051 ]. In our range of study, we hope to 
observe some traces of t his percolation as well as saturation properties observed in other moving oscillator 
systems at high speeds [Fujiwara et al . 2011bl ]. 



4 Prignano et al. 



It is also important to notice that we have kept the dynamical evolution of the units at its maximal 
simplicity since we are mostly interested in the interplay between motion (and hence construction of a 
dynamical network) and internal dynamics and how synchronization emerges as a collective property of 
the system. 



3. Results 

We present in figure [2] the results of our simulations. A preliminary observation points out that the roles 
and importances of V and r change throughout our map. 

For high enough values of r (r > r c ) , the synchronization time, T sync , is almost unaffected by the values 
of the speed V. Although this time is dependent on r, its range of possible values is much narrower (by 
orders of magnitude) than below the value, r c , that characterizes the static percolation transition. Actually, 
below this critical value the velocity plays the crucial role since for a fixed value of r the synchronization 
time changes by a factor of 10 4 . 
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Fig. 2. (Colors online) On the left: Heat map of the synchronization time as a function of r and V; in the top picture a large 
region of the parameters space has been considered, while the bottom figure is a zoom on the most critical region (r < 0.08, 
V < 0.25). On the right: Profiles of the efficiency of the system rT sync against r, for several values of V. 



To bring further insight to the map, we show a profile of the efficiency of the system rT sync dependence 
on r for various velocities. This efficiency variable balances the range of interactions (which has an "energy 
cost" related with the number of shots to different neighbors and their strength) and the number of shots 
that the system needs to synchronize (related with T sync ). 

It seems clear that mobility of agents helps to minimize energy consumption but as we approach the 
critical percolating value of r c , we observe that a range of behaviors emerge. On one hand, we observe for 
high velocities that the efficiency of the process remains roughly constant independently of r, since the 
extreme mobility of the agents compensate the reduced range of its interactions and successfully diffuses 
the sync hronization process around the sy s tem, a process that is equivalent to the observation in other 
settings (Frasca et all . 2008 : Fujiwara et all l2011b1 : lBaronchelli fc Diaz-Guileral . l201ll ]. O n the other hand, 
if the mobility of the agents is reduced, then the path to synchronization is more lengthly as well as more 
energy consuming. Synchronization is still possible (below the percolation static limit r c ) but the time to 



Synchronization of moving integrate and fire oscillators 5 



achieve it grows very fast, even resulting in an effectively infinite tim^]. 

These results lead us to identify different regimes and the consequent transitions between them. On 
the one hand, we find a "topological" transition at the critical static value r c , since above it there is 
basically no velocity dependence whereas below r c the influence of V is determinant. On the other hand, 
below r c , where synchronization is made possible by the mobility of the agents, we can identify a transition 
(depending on both r and V) that separates two dynamical regimes. 



4. Mechanisms 

Up to this point, we have studied the efficiency of our system and detected several regimes strongly 
dependent on two factors, the mobility of our agents and its range of interaction. 

In this part of the paper we want to study the distinct mechanisms that the system uses in its path to 
synchronization. To this end, we introduce an order parameter rj(T) = (cos (2ircj)(T))) that is an increasing 
function that measures the overall synchrony of our system, ranging from a uniform phase distribution of 
our oscillators (rj(T) = 0) to complete synchronizations (j](T) = 1). 

To couple our V and r control parameters we also introduce the number of distinct interactions per 
oscillator iV* (accumulated encounters with different agents by a oscillator i). This value is bounded (on 
average) between a minimal starting value of (iV c )o — {v — l)iV and a maximal value of (-^Vc)max = 

N-l and 

it provides information about the evolution of our system's synchronizing mechanism. Since the bounding 
of this value depends on r, we introduce a normalized magnitude x> 

(N c )(T)-(N c )q 

X{T) = (N c ) max - {N c)o W = *2>» < 4 > 

that is a quantification of the mixing of our system. When the mixing is minimal ((iV c ) = (N c )q) x 1S 0- 
As the system mixes, i.e. the oscillators increase its average number of contacted neighbors, it can grow 
up to its maximum value x = 1; i- e - (N c ) = {N c )max- 

We have calculated the time evolution of both r\ and % for different values of V € {0.1, 1, 100} and 
r G {0.05,0.1,0.2}. In figure [3] one can see the evolution of both parameters measured at the same time 
instants (and sufficiently averaged over enough realizations) together with the difference between them 
over time n{T) — x(T) that give insight about the topological evolution of our system as it synchronizes. 
The figures show three clearly distinct patterns. 

For high velocities we observe in fig. [3^) a gradual increase of the order parameter and a minor influence 
on r at fixed V, indicating that the synchronization emerges evenly on the system in a global fashion, due to 
the quickly changing topology of the network (neighbors of a given oscillator change rapidly). This regime 
(which we call diffusive) requires the inter-contact of the majority of the system, but this circumstance 
is rapidly achieved due to the strong mobility of the agents. In fact, this regime is optimal as far as 
the synchronization time is concerned, since th e interactions a r e more effective. These conclusions were 
obtained for populations of Kuramoto oscillators [Fujiwara et al . 2011b ]. for which this regime corresponds 
to the region of validity of the FSA. 

In the opposing case, when velocities are small enough, this behavior is completely lost and a step 
function appears, indicating that the slow mobility of the agents allows them to synchronize locally with 
their neighbors creating islands of synchrony. The sudden increase of the order parameter occurs at a 
regular pace, fact that points out that whenever the islands are disbanded (change of neighbors), they still 
transmit the local synchrony to the neighboring groups, mechanism that allows for system synchrony while 
keeping x i n small values. The initial height of the steps is dependent on r and decreasing as x(T) grows 
due to the limited range of rj(T) available states. 



^n our simulations some realizations of the experiment did not reach synchronization, even working with a reduced number 
of oscillators and the described boundary conditions, fact that induces us to think in this direction. 

2 Note that the average is calculated upon a firing event by the reference oscillator, hence our order parameter is an average 
of the phase difference of the other oscillators with respect to this one. 
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Fig. 3. (Colors online) Panel a): r\ against x f° r several values of r and V . Panel b): the difference between the two control 
parameters (77 — \) as a function of rT. Letters [D], [L] and [B] stand respectively for "diffusive", "local" and "bounded" 
regimes. The values of rj and x a t each time instant have been calculated averaging over 1000 realizations. 



Finally, as we decrease r approaching the critical value r c , we observe a transition from a local to 
a bounded regime, where the synchronizing time is so long that again allows for the interaction of the 
majority of the agents among themselves upon synchronization time (due to the bounded nature of the 
system). In this regime, the range of interaction is very reduced, and so is the size of the clusters, so an 
agreement between the multiple clusters created (if any) comes after almost all the system has interacted. 
Consequently, the increasing of rj with x is slower (many small steps, see fig. [3h.) ) while the final value x 
becomes larger (fig. E]). 

In fig. [3b) we provide an explicit time evolution of the difference rj(T) — x(T) in order to make the 
three regimes and the influence of r better identified. 

It is interesting to study the final mixing of our system upon synchronization as shown in figure [H 
This value x(T = T sync ) = 

Xsync together with T sync characterize the evolution of the system towards the 
synchronized final state. These features depend on both r and V . At fixed velocity the final mixing of the 
system decreases as the interaction range grows. This is caused by the fact that although an increased r 
induces more mixing (as oscillators find new neighbors more easily) it also drastically reduces (below the 
critical values r c ) the synchronizing time T sync thus reducing the chances of encounters between different 
oscillators. Above r c we observe a saturation of the values as the dependence of T sync in r and V is 
practically lost. 

Figure H]b) shows a change between the relation pattern of Xsync and efficiency rT sync ranging from the 
diffusive regime (mixing independent of rT sync ) to the curves where both the local and bounded regimes 
are shown. It is important to note the similarity of the observed shapes (for a wide range of V values) 
where we find the local phase concentrated around the minimum value of Xsync that gradually grows in a 
power law fashion as the performance of the system decreases (it consumes more energy to synchronize). 

The introduction of the new parameter Xsync allows us to present a phase diagram of our system 
relating the overall performance (in terms of efficiency) with the synchronizing mechanism used (figure [5]). 
We identify the diffusive regime in the zone of high velocities V ~ O(10) where both values of x an d rT sync 
are almost independent of r. This zone falters into the bounded one as V and specially r decrease, where 
both the efficiency and the mixing of the system is reduced. Finally for small enough velocities the local 
zone is clearly visible with low values of system mixing. From the map one clearly observes that the most 
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Fig. 4. (Colors online) On the left: the average value of \ at the synchronization time (xsync) as a function of r for several 
values of V. On the right: Xsync as a function of rT sync for several values of V. All the means have been performed over 200 
realizations. 



Xsync Xsync 




0.05 

Fig. 5. (Color online) On the left: The efficiency rT syn c as a function of r and V . The heat map of Xsync has been superimposed 
to the surface rT sync {r, V). On the right: the same quantities, zoomed in the most sensitive region of the parameter space. 

beneficial synchronizing mechanisms (in terms of energy consumption) are the diffusive and local ones. 
5. Conclusions 

Following the recent literature on complex systems, one of the hottest topics is the relation between 
dynamics and topology of interactions. In particular, there many evidences that the patterns of interaction 
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change rapidly with time, thereby completely altering and conditioning the dynamical properties of the 
system. 

Here we have proposed a framework in which agents move on a plane and are allowed to interact in 
a pulsing way. Each agent, representing a phase oscillator, moves at a common velocity and changes its 
internal phase with a common period. When this internal phase reaches a threshold value, the oscillator 
"fires", thus resetting its own phase as well as its orientation. At the same time, this oscillator interacts 
with the neighbors within a certain distance by changing their phases. Notably, this interaction setting 
makes the system reach a final synchronized state in which all oscillators fire at unison within the same 
fast scale. 

Keeping all geometrical parameters constant, we analyze the system behavior by changing solely the 
velocity of the agents and the interaction range. Notice that the motion of the agents is what ensures that 
the system will be able to synchronize. If the agents are static, a minimal fixed topology will be required 
to connect them (what is called a giant component in network terminology). In contrast, a population of 
moving agents (even when the interaction range is small) will eventually synchronize. 

We measure the time needed for the system to synchronize as a function of the two relevant agent 
parameters, the velocity and the fraction of population they interact with. Our model can be applied, 
for instance, to the field of wire-less communications by introducing the performance, which stands for 
the total number of signals emitted by the population to reach the synchronized state. Our simulations 
show that the time required for synchronization ranges widely, depending on the speed of the agents. As 
expected, the optimal performance is achieved when agents move very fast, irrespective of the interaction 
range; however, synchronization in systems with low mobilities dramatically depend on the interaction 
range. 

Finally, we have focused on the mechanisms that the system follows along its path to synchronization. 
We have introduced a new order parameter, that stands for the fraction of different units each oscillator has 
interacted with. This order parameter complements performance, and a novel phase diagram enables us to 
identify three different characterizable mechanisms leading to synchronization: i) the diffusive mechanism, 
where the system very quickly reaches synchronization by means of extremely effective interactions of 
each oscillator with a large fraction of the population; ii) the local mechanism, where agents only interact 
with a small population fraction, but local synchronization is sufficient to lead the system to the globally 
synchronized state; and finally iii) the bounded mechanism, characterized by very slow motion and short 
range of interaction, which allow a high degree of accumulated mixing during a very long synchronizing 
time. 

The identification of these mechanisms that relate mobility and interaction to synchronization will 
undoubtedly be crucial in similar models of populations of moving agents. Indeed, understanding these 
phenomena will help to design optimal protocols to dissect more realistic settings, as well as to predict 
their behavior, by defining their interaction rules. 
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